Local vortex nucleation and the surface mode spectrum of large condensates 
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A combination of analytical and numerical approaches ob- 
tains the complete dispersion curve for surface excitations in 
a condensate held in a plane linear potential. This improve- 
ment on previous approximate results yields an accurate for- 
mula for the local Landau critical velocity for vortex nucle- 
ation at the surface of a sufficiently large condensate, which 
agrees very well with recent experiments. The dispersion 
curve for surface modes at a hard wall potential is also pre- 
sented, for contrast. 
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Critical velocities for vortex nucleation in dilute 
trapped Bose-Einstein condensates have recently been in- 
vestigated theoretically and experimentally dj-^. 
It is clear topologically and energetically that unsta- 
ble surface excitations are the seeds from which vor- 
tices grow, and this fixes the characteristic scales of 
the problem. Quantitative comparison between theory 
and experiment, however, requires more than order-of- 
magnitude estimates. A precise but general theory may 
be based on the idea that the Thomas-Fermi (TF) sur- 
face of a large condensate enters a 'bulk regime' in which 
the physics of vortex nucleation is local. A merit of this 
local theory is that it can be applied to generic conden- 
sate flows, which cannot naturally be described in terms 
of global surface modes of definite angular momentum. 

When vortices are produced 'gently' enough that one 
can consider them to grow from a quasi-steady state, the 
local theory of vortex nucleation proceeds in three steps. 
The first is determining the condensate density and ve- 
locity field in the vortex-free steady state, as a function 
of experimental control parameters such as atom num- 
ber and stirring frequency. The second is determing the 
local critical velocity over the TF surface (which in the 
presence of stirring beams may be multiply connected), 
and so assessing when and where the surface velocity 
field may be locally supercritical. The third is com- 
puting the rate at which supercritical surface excitations 
grow into nonlinear vortices. The first step is in general 
a difficult problem in three-dimensional hydrodynamics, 
though in some cases it can be performed analytically. 
The second step is the subject of this paper. The third 
step again depends on experimental details, such as ki- 
netics of the thermal cloud, strength of stirring potential, 
and timescale over which control parameters are changed. 



Whatever issues arise in steps one and three, the Landau 
criterion | pO{ is universal in step two: an energetically un- 
stable surface mode must appear for vortices to enter the 
cloud. For sufficiently large condensates, such as some 
in recent experiments, the local Landau critical velocity 
may be determined by numerically solving one universal 
equation. 

Near the Thomas-Fermi (TF) surface of a large 
trapped condensate, the trapping potential is approxi- 
mately linear. If the condensate is large enough, and 
does not have too extreme an aspect ratio, local physics 
should also be insensitive to the curvature of the TF sur- 
face. So, following Al Khawaja, Pethick and Smith |Tl| ] 
(hereafter AKPS), one can approximate a generic TF 
surface by replacing the actual trapping potential with 
a linear ramp V = Fx for constant F. The equations 
simplify by rescaling from physical to dimensionless vari- 
ables, defining x = S^'^Xph^ t = r^^tph, ip = dVSTratljph, 
where S and t arc the characteristic surface length and 
time scales 
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M is the particle mass, a is the s-wave scattering length, 
and ipph is the macroscopic wavefunction. (The chemi- 
cal potential factor e"'''*^'' has been extracted from ipph, 
which is normalized so that iV'phl^ is the density of par- 
ticles in the condensate.) In the dimensionless variables 
the Gross-Pitaevskii equation then reads 
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where the origin of x is at the TF surface. Setting 
dttp = determines the stationary background solution 
ip = f{x), where / is taken to be real. It may be found 
to good accuracy numerically . 

The dispersion relation for surface excitations may be 
found by solving the Bogoliubov equations for small per- 
turbations ip f + rj 5ip, with infinitesimal rj and 

= u{x) e'('J''-"*) V* {x) e-'(9''-"*) . (3) 

Here the excitation's wave vector is taken to lie along 
the y axis, and the dimensionless frequency U, = lot and 
wave number q — kS are defined. Note that far outside 
the TF surface, for x — s- oo, f{x) cx e~i^ — > 0, and 
hence 



1 



lim Si: (X e'('?y-"*) exp[--x^/2 + (f^ - q^)x^/^ 

X — *oo 3 



lim Tp (X e 3^ 



Since it turns out that fl > q^, this means that the growth 
of 77 does not just lead to vortex formation, but actually is 
the approach of vortices to the TF surface from infinity. 
It is convenient to solve for the dispersion relation 
by iterating the Bogoliubov equations to obtain the 
fourth order eigenvalue problem for s{x) = u — v: 
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Discretizing x and representing the second derivatives 
with finite differences converts this fourth order ODE di- 
rectly into a matrix diagonalization problem, which can 
be solved numerically. Taking ^(q) to be the square 
root of the lowest eigenvalue fl^ for given q, one finds 
the function shown in Figure 1. (For any q there is 
obviously a discrete spectrum of fl^, but the critical ve- 
locity is set by the lowest branch.) Numerical solution 
becomes difficult at small q, where even the lowest lying 
eigenstate s{x) extends far enough into the condensate 
that a large grid is needed to contain it; but fortunately 
the problem remains easily tractable until the numeri- 
cal curve has already converged onto the ^/2q asymp- 
totic behaviour as g ^ 0, found analytically by AKPS 
| pj| . At large q, simple perturbation theory implies that 
n{q) ^ q^ + Eg2 + 0{q'^^), where Eg2 = 1.17 is the lowest 
eigenvalue of i?2. (This is the same as the 'Popov ap- 
proximation' of neglecting the anomalous average in the 
Bogoliubov equations for u and v.) Figure 1 shows that 
the exact dispersion curve approaches these asymptotic 
forms closely for g < 0.3 and for q > 1.5. 

For comparison. Figure 1 also shows the approximation 
of AKPS UM- which in our dimensionless units is 
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and was derived to give the first corrections at small q 
beyond the leading term il^ ~ 2q. According to Figure 
1, it clearly does this very well, tracking the exact disper- 
sion curve closely for q < 0.6. This approximation was 
never intended to be accurate at larger q, however, and it 
obviously becomes pathological for g > 1; its application 
to vortex nucleation in Q is therefore only qualitatively 
accurate. The recent on the other hand, is based on 
the high q regime, in the sense that it fits numerical data 
to the functional form q^+ const., and so underestimates 
Vc by about 15%, because this form is not really valid near 
qc- For comparison with (^), in the range < g < 2 the 
numerical result shown in Figure 1 is indistinguishable, 
on the scale of the plot, from 



This formula becomes inaccurate at 5 ^ 2, where it is 
better to use the asymptotic Popov result. 
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Figure 1. Dimensionless surface mode dispersion relation 
(solid curve). This curve is numerical for q > 0.1; for 
g < 0.1 it is simply ^/2q. The tangent to fl{q) through 
the origin has slope 2. The approximation of Al Khawaja, 
Pethick and Smith is shown for comparison. Dotted curves 
show the asymptotic behaviours ^/2q and + Eg2 (see text). 

The main result which we can derive from our numeri- 
cal knowledge of ^{q) is an accurate value for the Landau 
critical velocity due to surface excitations, which in units 
of 5/t is simply the slope of the tangent to ^{q) through 
the origin. Numerically this slope is indistinguishable 
from 2. This leads to the expression in physical units 
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This formula can immediately be obtained as an order 
of magnitude estimate, of course; but the fact that it is 
precisely accurate is nontrivial. Despite its simple ap- 
pearance, I have been unable to derive this result analyt- 
ically, and it may not actually be exact; but it is accurate 
at least to within a few tenths of one percent. Another 
result which we obtain is the wavenumber qc at which 
fl{q)/q is minimized, which is approximately qc = 0.89. 
This number is close to 2~^/^; but since ^{q) has rather 
small curvature near q ~ qc, the numerical uncertainty 
in qc is too large to make this identification more than a 
conjecture. 

The applicability of the local surface critical velocity 
(P) to recent experiments may be assessed as follows. 
Since 5 is the characteristic depth to which surface modes 
of q ~ 1 penetrate the condensate, the bulk regime re- 
quires 5/R = e <^ 1 (where R is the TF radius of the con- 
densate at its rotational equator) so that surface modes 
are indeed confined to within the range over which the 
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harmonic potential is accurately linear. Since ttS is also 
the characteristic surface wavelength of the critical ex- 
citations, neglecting the curvature of the TF surface re- 
quires that S is much smaller than either radius of cur- 
vature. In prolate traps rotated about the long axis, the 
shortest radius of curvature is just R; but in oblate traps 
the curvature in the axial (z) direction is greater, and we 
need the more stringent condition e = {lo16) / {ui'^ R) <^ 1. 

For a rotationally symmetric harmonic trap, the quan- 
tity {qcR/S — 1) may be compared to the multipolarity 
Ic of the rotationally symmetric surface mode, and the 
surface critical velocity can be converted into the dimen- 
sionless critical rotation frequency as a fraction of radial 
trapping frequency, given in the Thomas-Fermi limit of 
large particle number N by 
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For experiments reported to date, computing F using 
the TF density profile for the reported numbers of atoms 
yields 



Experiment 


A^ 


£ 


qcRI5 


X 


Xobs 


MIT stiff |7] 


1.5 X 10^ 


0.045 


19.6 


0.30 


0.29 


MIT weak || 


5 X 10^ 


0.030 


29.8 


0.24 


0.27 


ENS [ 


6| 


3 X 10^ 


0.11 


8.4 


0.46 


0.65 


JILA 


1 


6 X 10'^ 


0.12 


30.7 


0.24 


0.3(5 ±3) 


Oxford M 


2 X 10"* 


0.72 


9.9 


0.42 


0.56-1- 



The values of e (which include oblateness factors of 4 
and 8 for the JILA and Oxford traps respectively) in- 
dicate that the local surface theory should describe the 
MIT experiments very well, and should apply to the JILA 
and ENS experiments with fair accuracy, but is not ex- 
pected to fit the measurements at Oxford (which were 
made at various eccentricities, not reflected in the table, 
and to which the global hydrodynamic analysis of 1^,?] 
is relevant instead). 

The local theory actually agrees even better with the 
MIT experiments than the table indicates, because with 
the stirring beam applied the MIT traps were signifi- 
cantly asymmetric, and for the weak configuration even 
anharmonic. Hydrodynamics in the weaker trap is there- 
fore not analytically tractable, but the stiff MIT trap re- 
mains harmonic with the stirring beam applied, and so 
the hydrodynamic solution for a rotating harmonically 
trapped condensate yields the non-perturbative den- 
sity and velocity fields, in the co-rotating frame. One 
can use this solution straightforwardly, for a given trap 
eccentricity and rotation frequency, to compute the local 
condensate velocity v{9, </>) over the whole TF surface, as 
well as the local surface force F(9, </>) (taking centrifugal 
effects into account), and hence the local critical velocity 
Vc{0,(t)). The results show that the local flow velocity 
first reaches the local critical velocity at the poles of the 



cloud's shortest axis, so that vortices will first form in 
these surface regions. One also finds that local critical- 
ity is reached in the stiff MIT trap at a slightly lower x of 
0.285. (Although angular momentum is not conserved in 
the asymmetric trap, the angle subtended by 2t:5 /qc at 
the vortex nucleation radius (smallest semi-axis) yields 
the effective I = 18 cited in |^.) 

The excellent agreement with the MIT experiments is 
likely due not only to the lower values of e, but also to 
the facts that the MIT experiments were most sensitive 
in detecting vortices (due to condensate size and a 'slic- 
ing' technique to improve visibility), and applied strong 
stirring fields (so that the timescales for vortex nucle- 
ation were probably shorter for a given degree of super- 
criticality). As suggested in 1^, the nucleation time in 
the weakly anisotropic ENS trap may be too long for 
visible vortices to be created, except in the large (and 
unstable j|]) velocity fields generated by rotations near 
the quadrupole resonance. Similarly, for the JILA sys- 
tem the X — 0-24 given by ^ is significantly lower than 
the lowest value 0.32 at which vortices may have been 
detected [Q. This discrepancy is enough to motivate 
further development of the local critical velocity theory, 
to determine whether ©(e) corrections due to surface cur- 
vature may have large constant factors. It is also worth 
noting, however, the significant experimental uncertainty 
due to the indirect detection of small numbers of vor- 
tices through changes in the condensate aspect ratio 
And since vortex nucleation in the JILA experiments is 
through equilibration with the rotating thermal cloud, 
one can estimate the rate of surface instability growth 
using quantum kinetic theory, and show that vortex nu- 
cleation is slow until significantly supercritical velocities 
are reached. (When the thermal cloud's mean velocity 
past the TF surface just exceeds Vc |15|, the highest gain 
over loss is for the mode with q ^ qc, and is propor- 
tional to (1 — exp[— 2/3(7c(f — fc)]), where (3 is the thermal 
cloud's inverse temperature in our dimensionless surface 
units. Since experimental temperatures of T ^ lOOnK 
are much higher than the temperature corresponding to 
the surface critical velocity, this again indicates that the 
timescale for vortices to form may be quite long until 
the critical rotation frequency has been significantly ex- 
ceeded.) Further study of such timescales is needed, but 
it appears that the JILA data do not demand any vortex 
nucleation mechanism other than surface mode excita- 
tion by the 'wind' of the rotating thermal cloud. 

More recent observations at MIT ||l^ using a narrow, 
stick-like stirring beam have obtained x ^-s low as 0.08, 
where Eqn. (^ predicts 0.24. As in the rotating stiff 
trap, however, the local critical velocity must be com- 
pared to the actual local fluid velocity created over the 
entire TF surface, including the surface at the stirring 
beam. Detailed analysis of the 3D flow induced by the 
Gaussian beam will require substantial numerical effort, 
and is beyond the scope of this paper; but simplified sim- 
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ulations show that fluid velocities near Gaussian beams 
can easily be many times the velocities of the beams 
themselves. (The well-known hydrodynamic speed-up 
factor of 2 for a hard cylindrical stirrer is enhanced be- 
cause the fluid density drops as it penetrates the beam.) 
Of course if the beams are too narrow in relation to their 
strength, the ratio of surface depth to curvature radius 
on the TF surface may not be sufficiently small for the 
zeroth order theory presented here to be accurate. On 
the other hand, if the S on the beam surface is actu- 
ally smaller than the healing length ^ = (47ra|V'p/iP)~^^^ 
near the beam, the linearized theory will also fail, as 
the beam potential approaches the limit of a hard wall. 
The surface mode spectrum of a hard wall can also be 
computed, if we can again neglect curvature. In this 
case we have an analytic solution for the stationary back- 
ground ip oc tanhx/^, and a multiple scale analysis valid 
at long wavelength yields lowest Bogoliubov frequencies 
uj^ EE (2M(^/h)uj = 2ke, + 0[{kS)% As shown in Fig. 
2, the result is that the surface mode dispersion curve 
uj(^(k^) for a hard wall lies below the bulk dispersion 
curve, but not by enough to lower the critical velocity 
below the local speed of bulk sound. 




Figure 2. Hard wall surface mode dispersion relation tj^(fc^) 
(solid curve). Results are numerical only for q > 0.1; be- 
low this point the solid curve is simply 2fc^. For all fc, the 
surface mode dispersion curve lies between 2fc^ and the bulk 
dispersion curve fc^y4 + fc2^. 

To summarize, in condensates of at least moderate size 
vortex nucleation at TF surfaces is governed by the lo- 
cal Landau critical velocity obtained from the dispersion 
relation for Bogoliubov surface modes. This velocity is 
h/[M5) to high accuracy (not just of that order). Ex- 
periments at MIT support this theory well. Failure to 
see vortices at such low velocities elsewhere is likely due 
to long nucleation times. For very abrupt potentials, 
the linear theory breaks down and the local surface mode 
critical velocity approaches the local speed of bulk sound. 
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The surface mode x = 0-40 cited in [g] is apparently 
based on the sum rule upper bound of F. Dalfovo and 
S. Stringari, Phys. Rev. A 63, 011601(R) (2001), which 
is not expected to be accurate at the high I = 30 corre- 
sponding to q — Qc for this system. 

Imposing conservation of energy and momentum in colli- 
sions that generate surface excitations, one obtains that 
a surface mode of dimensionless wavenumber q can only 
be excited by incident thermal atoms whose momen- 
tum component in the direction of the surface mode 
is qi = [il{q) + q'^]/{2q). The same surface mode can 
only be de-excited by incident thermal atoms of qi — 
[f2(g) — g^]/(2q). For a Boltzmann distribution of ther- 
mal atoms P oc exp[—/3{qi — v/2)^], therefore, the ratio 
between rates of gain and loss can only exceed unity for 
any surface mode q if v > Vc. So even though the thresh- 
old velocity for an incident atom to excite surface modes 
at qc is higher than Vc, due to atomic recoil, the critical 
velocity for a thermal cloud of non-interacting atoms is 
still exactly v^. 
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